Cite this Article
Wang Yu, Zhang Rongpei, Wang Zhen, Han Zijian. Turing pattern in the fractional Gierer–Meinhardt model. Chinese Physics B, 2019, 28(5): 050503  
Permissions
Turing pattern in the fractional Gierer–Meinhardt model
Wang Yu1, Zhang Rongpei1, †, Wang Zhen2, Han Zijian1
School of Mathematics and Systematic Sciences, Shenyang Normal University, Shenyang 110034, China
College of Mathematics and Systems Science, Shandong University of Science and Technology, Qingdao 266590, China

 

† Corresponding author. E-mail: rongpeizhang@163.com

Abstract
Abstract

It is well-known that reaction–diffusion systems are used to describe the pattern formation models. In this paper, we will investigate the pattern formation generated by the fractional reaction–diffusion systems. We first explore the mathematical mechanism of the pattern by applying the linear stability analysis for the fractional Gierer–Meinhardt system. Then, an efficient high-precision numerical scheme is used in the numerical simulation. The proposed method is based on an exponential time differencing Runge–Kutta method in temporal direction and a Fourier spectral method in spatial direction. This method has the advantages of high precision, better stability, and less storage. Numerical simulations show that the system control parameters and fractional order exponent have decisive influence on the generation of patterns. Our numerical results verify our theoretical results.

PACS: 05.30.Pr;05.65.+b;82.40.Ck;02.70.Hm
Keyword:Turing patterns;fractional Gierer-Meinhardt model;Fourier spectral method
1. Introduction

Reaction–diffusion systems have been used to describe pattern formation. For example, Turing found that diffusion could cause instability of a spatially uniform state and lead to the formation of spatial patterns.[1] This diffusion-driven instability or Turing instability has attracted attention of scientists in many areas, such as biology,[2] physics,[3] neuroscience,[4] optics,[5] chemistry,[6] and geology.[7] Typical examples of pattern formation are found in the Schnakenberg model,[8] the chloride-iodide-malonic acid reactive model,[9] the Gray–Scott model,[10] and the Gierer–Meinhardt model.[11] These mathematical models have the general form

where u demonstrates the autocatalytic behavior and is called activator, while v hinders the activator autocatalysis and is called inhibitor. Du, Dv are constant diffusion coefficients and f, g describe the reaction kinetics.

The normal diffusion in system (1) corresponds to Brownian motion which has finite moments in both waiting time distribution and jump size distribution. However, particles do not always carry out the recent jumps but can instead wait and perform long jumps between successive jumps. When the finite moments conditions in the waiting time distributions and jump size are not satisfied, anomalous diffusion of particles occurs.[12,13] There are two types of anomalous diffusions: superdiffusion and subdiffusion. Superdiffusion is faster than the normal diffusion, and the corresponding jump size distribution has infinite moments, or alternatively, Lévy flights.[14] Subdiffusion is slower than the normal diffusion, and the corresponding waiting time distribution has infinite moments. Both types of anomalous diffusions play important roles in biological, physical, and chemical processes. It has also been shown that anomalous diffusion[15] in one-dimensional systems leads to the occurrence of anomalous heat conduction. The spiral wave and chemical turbulence were observed in the nonlinear dynamics of the oscillating reaction diffusion patterns in Ref. [16]. In Ref. [17], the pattern formation in an Oregonator model with superdiffusion was studied.

The pattern formation occurring in anomalous diffusion cases is depicted by space fractional reaction–diffusion equations

where , depict the fractional superdiffusion. The initial conditions are chosen as small perturbations away from the steady state of both substances. Zero-flux boundary conditions are chosen in the form , because patterns are evolved only by self-organization. To examine the diffusion-driven instability conditions for the reaction–diffusion equations (2), we need to find the eigenvalue problem and evaluate its characteristic polynomial.

There is no universal and effective way to obtain the exact solutions for most fractional differential equations. Consequently, numerical solutions for the fractional differential equations are needed. At present, the numerical approaches and supporting analysis of fractional-order differential equations are a fruitful research topic. Li and Xu[18] considered a spectral method for the weak solution of the space–time fractional diffusion equation. Khader et al.[19,20] chose Chebyshev and Legendre Galerkin methods to discrete the fractional advection-dispersion equations where the spatial derivative was selected as the Caputo derivative. Hanert[21] proposed the Chebyshev spectral element method to solve the fractional Riemann–Louiville advection-diffusion equation. Recently, Zayernouri and Karniadakis[22] presented a collocation method based on fractional Lagrange interpolants for solving the steady-state and time-dependent fractional partial differential equations.

In this paper, we use the Fourier spectral method as an efficient alternative approach in spatial discretization. In terms of time discretization, we use the explicit fourth-order time stepping method, exponential time differencing Runge–Kutta method, which improves the stability and computation efficiency. The exponential time differencing method was proposed by Cox-Matthews in Ref. [23] and improved by Kassam-Trefethen in Ref. [24]. Specifically, the main contributions of this work are as follows:

(i) We derive the conditions of Turing instability for the fractional Gierer Meinhardt system by using linear stability analysis. The various numerical examples verify the linear stability analysis results which show that the system control parameters and fractional order exponent have decisive influence on the generation of patterns.

(ii) We apply Fourier spectral methods to fractional-in-space reaction-diffusion equations. The main advantage of this method is that it gives the full diagonal matrix expression of fractional operators, which has spectral convergence and does not have to take the fractional power into consideration.

(iii) We use the exponential time differencing Runge–Kutta method, which has a much better performance in computer time.

The rest of this paper is organized as follows. Section 2 gives the spectral decomposition of fractional Laplacian operator, then the mathematical mechanism of the fractional pattern and the conditions inducing the Turing instability are determined. In Section 3, we apply the Fourier spectral method in spatial discretization and exponential time differencing scheme in time discretization for solving the reaction diffusion systems. In Section 4, the numerical examples of fractional Gierer–Meinhardt model are given to verify the result of our analysis. Finally, we present a summary of this paper.

2. The mathematical mechanism of Turing pattern

In this section, we first give the definition of the spectral decomposition of fractional Laplacian operator. Then, we perform the linear stability analysis for the general reaction–diffusion system (2). Finally, a set of conditions for Turing instability are derived for the fractional Gierer–Meinhardt model.

Spectral decomposition plays a crucial role in the definition of the fractional Laplacian operator, , . We first define the spectral decomposition for Laplacian operator . Suppose has a complete set of orthonormal eigenfunctions satisfying homogeneous Neumann boundary condition on a bounded region , with corresponding eigenvalues i.e., on . We can define the space as follows:

For any , the fractional Laplacian operator can be defined as

For the zero flux boundary condition, it can be easily obtained that there exist , in one dimensional case , and , in two dimensional case with j, It can be observed that has the same interpretation as in terms of its spectral decomposition.

Assume that the homogeneous steady state corresponds to the solution of , Omitting the diffusion terms in system (2), the spatially homogeneous system becomes

Take linearization for the nonlinear ordinary differential equations (5) around the steady state and let , we obtain
The homogeneous steady state (U,V)=(0,0) is stable when the real part of the eigenvalues of Jacobian matrix is negative, i.e., .

Take linearization for the nonlinear fractional reaction diffusion equations (2) around the steady state and let , we obtain

To examine the diffusion-driven instability conditions for the fractional reaction–diffusion systems, we need to compute the eigenvalue for the following equation:

with homogeneous Neumann boundary condition.

Assume that U, V are decomposed by the complete set of orthonormal eigenfunctions and . Substituting it into Eq. (7) leads to the following linear algebraic systems:

where for 1D case and for 2D case. The necessary and sufficient condition of non-zero solutions for the homogeneous linear systems (8) is

Now we apply the linear stability analysis in a specific activator–inhibitor system, Gierer–Meinhardt system. In 1972, Gierer and Meinhardt proposed a prototypical activator and inhibitor system to model the generating mechanisms of biological patterns.[25] The Gierer–Meinhardt model is widely used in the study of some basic phenomena in the process of morphogenesis.[26] To our knowledge, there are few published works on the two-dimensional spatial fractional Gierer–Meinhardt models. In the present paper, we consider the following fractional Gierer–Meinhardt model:

where the parameters are taken as and η=0.1. It can be easily obtained that the stead state for fractional Gierer–Meinhardt model (9) is .

The Jacobian matrix is , and the characteristic polynomial for the fractional Gierer–Meinhardt model is

The corresponding eigenvalues are

If the Turing pattern is formed as a result of system evolution, there exists a range of wave numbers μ such that . We call the function the dispersion curve, in which the information about existence of nonuniform stationary patterns can be obtained.

For different κ, the different dispersion curves are plotted in Fig. 1. For every κ, four dispersion relations are yielded for different combinations of fractional order α1 and α2. In our computation, we find that the eigenvalue λ is real. Figure 1(a) shows the dispersion curves for κ=0.0162 with fixed α1=2 and different order α2. It can be observed that changes from positive to negative with decreasing α2, which means that the homogeneous state finally becomes stable as time evolves for α2=1.7.

Fig. 1. The relation of eigenvalue λ with μ for (a) κ=0.0162, (b) κ=0.0128, and (c) κ=0.008.

Figure 1(b) shows the dispersion curves for κ=0.0128 with different values of α1 and α2. Under the circumstance of α1=2 and α2=1.8, and the system is stable. While in other three cases, the stationary pattern is formed as system evolves because of .

Fix the fractional order α2=2 and and take α1=2, 1.9, 1.8, and 1.7, the dispersion curves for κ=0.008 are plotted in Fig. 1(c). One can notice that the Turing instability occurs as α1 decreases. Thus, we can conclude that the fractional orders α1 and α2 have significant effects on the dispersion relations. Turing patterns can be controlled by changing the value of the fractional Laplacian in Eq. (2).

3. Numerical method
3.1. Space discretisation

The two-dimensional computation domain is discretized into tensor product grids . The spatial mesh sizes are chosen as and . We consider a finite number of orthonormal trigonometric eigenfunctions

where and .

The solution u at can be expanded by the truncated series

The coefficients can be efficiently computed by fast direct discrete cosine transforms (DCT)
By applying the discrete cosine transform to the fractional reaction diffusion equation (2), we obtain two-dimensional representation in spectral space

The linear term of the ordinary differential equations is diagonal. Moreover the nonlinear term can be evaluated in physical space and then transformed to spectral space. Here, and are the double cosine transforms of the chemical concentrations u and v, respectively. and are the double cosine transform coefficients of the matrix and , .

3.2. Time discretization

We consider the time discretization for Eq. (12). This method is similar to the integral factor method. By multiplying Eq. (12) by the integration factors and , we obtain

Here we take the first equation in Eq. (13) as an example to apply the exponential time differencing formula of Runge–Kutta type. The process for the second equation is similar. We integrate the first equation in Eq. (13) over a single time step of length τ (from to ), getting
where the matrix is formed by and elements in are . Many existing exponential time differencing methods are used to approximate the integral terms in this formula. Matthews utilized exponential time differencing methods based on Runge–Kutta time step in Ref. [23], which is called exponential time differencing Runge–Kutta method. Set , then the fourth order exponential time differencing Runge–Kutta scheme is given as follows:
where the stages , i=1,2,3 are given as
and the functions , i=1,2,3 are defined as
The final solution can then be obtained by applying the cosine inverse transforms for .

4. Numerical experiments

Example 1 We consider the following fractional reaction–diffusion equation on :

subject to the homogeneous Neumann boundary condition. The initial data is given as and right-hand side . For the convenience of comparison, we solve this problem by the Fourier method and finite difference method (FDM) proposed in Ref. [27]. We demonstrate the efficiency and accuracy of the two methods by reporting the -norm error
where u is the numerical solution and is the reference solutions calculated by evaluating the fractional equation (17) with 212 Fourier modes. The convergence results in the -norm are presented in Table 1 at t = 1 with D = 1 and α=1.8. As seen in the table, the Fourier method is more accurate than the FDM. The CPU time is also give in Table 1 for comparison between the methods. It can be concluded that our method has a much better performance in execution time.

Table 1.

The error and CPU time for Fourier and FDM methods in Example 1 at T = 1 with τ=0.001.

.

Example 2 Consider the fractional Gierer–Meinhardt model (9) on with homogeneous Neumann boundary condition. The parameters are taken as and η=0.1. We consider three cases of κ=0.0162, κ=0.0128, and κ=0.008. The initial conditions are chosen as

In the following numerical tests, the spatial mesh sizes are chosen as . The time step is chosen as .

We first consider the case of κ=0.0162 and fix α1=2. From the discussion on the stability of the Gierer–Meinhardt system in Section 2, Turing instability will occur for α2=2, 1.9, and 1.8. And the systems will become stable for α2=1.7. We plot the solution of u at different time in Fig. 2 for different α2. As seen in Fig. 2(a), the solution u approaches the steady state for α2=1.7. The steady state breaks up and various patterns can be observed by increasing α2. The stripe pattern can be observed in Fig. 2(b) for α2=1.8. If we keep increasing α2, then a chain cluster of Turing spotted patterns initially evolve with α2=1.9 and 2; as seen in Figs. 2(c) and 2(d). As the simulation time increases to 700, we obtain pure Turing spot patterns. The stationary state spatial structures have a period which grows as α2 increases. These results confirm the presence of linear stability analysis for the fractional Gierer–Meinhardt model; that is, Turing instability occurs when the real part of eigenvalue is positive.

Fig. 2. Turing pattern for κ=0.0162, α>1=2, (a) α>2=1.7, (b) α>2=1.8, (c) α>2=1.9, and (d) α>2=2.

Next we consider the case of κ=0.0128. Figure 3 shows the scenarios of pattern formation with different diffusion indices. As the dispersion curve shows in Section 2, the systems will become stable for α1=2 and α2=1.8, see Fig. 3(a). While Turing patterns appear when the parameters are set as α1=2, α2=1.9, α1=2, α2=2, and α1=1.9, α2=2. With increasing α2 from 1.9 to 2, a stripe pattern turns into a spot pattern. Fix α2=2 and decrease α1, we can observe that the period of spot becomes larger.

Fig. 3. Turing pattern for κ=0.0128, (a) α>1=2, α>2=1.8; (b) α>1=2, α>2=1.9; (c) α>1=2, α>2=2; and (d) α>1=1.9, α>2=2.

Finally, we consider the case of κ=0.008. Figure 4 shows the scenarios of pattern formation with . The Gierer–Meinhardt model becomes stable as time advances when α1=2 and α2=2. The component of u is shown in Fig. 4(a). When α1=2,1.9, 1.8 and α2=2, Turing patterns then appear. The common feature of Figs. 4(b)–4(d) is that the patterns are formed at the initial stage of diffusion, and remain unchanged. The spots become sparse with decreasing α1.

Fig. 4. Turing pattern for κ=0.008: (a) α>1=2, α>2=2; (b) α>1=1.9, α>2=2; (c) α>1=1.8, α>2=2; and (d) α>1=1.7, α>2=2.

Remark 1 In this paper, we use the same initial conditions (18) as integer case[26] to generate the Turing pattern. The other initial conditions can generate similar results. We also compute the fractional Gierer–Meinhardt model with the random initial conditions. The stable pattern with κ=0.0162, α1=2, and α2=1.7, 1.8, and 1.9 are shown in Fig. 5. Compared to Fig. 2, it can be found that the solution u approaches the steady state for α2=1.7. Moreover, the stripe pattern also preserves for α2=1.8 and α2=1.9, leading to the spotted patterns.

Fig. 5. Turing pattern for κ=0.0162 with random initial condition: (a) α>2=1.7, (b) α>2=1.8, and (c) α>2=1.9.
5. Conclusion

We investigated the mathematical mechanism of the fractional Turing pattern by linear stability analysis. Our analyses help to explain the generation of the Turing pattern and calculate the threshold of the instability. We combined the Fourier spectral method for space discretization and exponential time differencing Runge–Kutta method for time discretization to simulate the fractional reaction–diffusion equations. Numerical experiments show that the fractional reaction–diffusion equations can exhibit the various characteristics of the Turing pattern. Besides the system control parameters, the fractional order exponent also decides the generation of Turing pattern.

Reference
[1] Turing A M 1952 Philos. Trans. Royal Soc. 237 37
[2] Baurmann M Gross T Feudel U 2007 J. Theor. Biol. 245 220
[3] Ouyang Q Swinney H L 1991 Nat. 352 610
[4] Liang P 1995 Phys. Rev. Lett. 75 1863
[5] Staliunas K Sánchez-Morcillo V J 2000 Opt. Commun. 177 389
[6] Epstein I R 1984 Nat. 307 692
[7] L’Heureux I 2013 Philos Trans. 371 20120356
[8] Iron D Wei J Winter M 2004 J. Math. Biol. 49 358
[9] Li Q Zheng C G Wang N C Shi B C 2001 J. Sci. Comput. 16 121
[10] Wang T Song F Wang H Karniadakis G 2019 Comput. Meth. Appl. Mech. Engrg. 347 1030
[11] Yang R Song Y 2016 Nonlinear Anal. Real World Appl. 31 356
[12] Gambino G Lombardo M C Sammartino M Sciacca V 2013 Phys. Rev. 88 042925
[13] Golovin A A Matkowsky B J Volpert V A 2008 SIAM J. Appl. Math. 69 251
[14] Metzler R Klafter J 2004 J. Phys. A: Math. Gen. 37 R161
[15] Li B Wang J 2003 Phys. Rev. Lett. 91 044301
[16] Nec Y Nepomnyashchy A A Golovin A A 2008 Physica 237 3237
[17] Feng F Yan J He Y F Liu F C 2016 Chin. Phys. 25 104702
[18] Li X Xu C 2010 Commun. Comput. Phys. 8 1016
[19] Khader M M 2011 Commun. Nonlinear. Sci. Numer. Simul. 16 2535
[20] Khader M M Sweilam N H 2014 Comput. Appl. Math. 33 739
[21] Hanert E 2010 Environ Fluid Mech. 10 7
[22] Zayernouri M Karniadakis G E 2014 SIAM J. Sci. Comput. 36 A40
[23] Cox S M Matthews P C 2002 J. Comput. Phys. 176 430
[24] Kassam A K Trefethen L N 2005 Soc. Ind. Appl. Mech. 26 1214
[25] Gierer A Meinhardt H 1972 Kybernetik 12 30
[26] Zhang R P Wang Z Wang Y Han Z J 2018 Acta Phys. Sin. 67 050503 (in Chinese) http://dx.doi.org/10.7498/aps.20171791
[27] Bueno-Orovio A Kay D Burrage K 2014 BIT Numer. Math. 54 937